!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine roe (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                           &
   ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, gm)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx
   real(8), intent(in) :: gm
   real(8), dimension(nx), intent(in) :: ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l
   real(8), dimension(nx), intent(in) :: ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r
   real(8), dimension(nx), intent(inout) :: fro, fmx, fmy, fmz, fbx, fby, fbz, fen

   integer(4) :: i, ip1, i_dim
   real(8), dimension(nx) :: ro_a, mx_a, my_a, mz_a, bx_a, by_a, bz_a, en_a
   real(8), dimension(8) :: var_c
   real(8), dimension(8, nx) :: eigenvalue
   real(8), dimension(nx) :: fro_l, fmx_l, fmy_l, fmz_l, fbx_l, fby_l, fbz_l, fen_l
   real(8), dimension(nx) :: fro_r, fmx_r, fmy_r, fmz_r, fbx_r, fby_r, fbz_r, fen_r
   real(8), dimension(8, nx) :: var_d
   real(8), dimension(8, 8, nx) :: eigenvector_l, eigenvector_r
   real(8) :: gmm1, vx, vy, vz, b2, v2, pr

!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0

   do i = 1, nx - 1

      ip1 = i + 1

      vx = mx_l(i) / ro_l(i)
      vy = my_l(i) / ro_l(i)
      vz = mz_l(i) / ro_l(i)
      b2 = bx_l(i) * bx_l(i) + by_l(i) * by_l(i) + bz_l(i) * bz_l(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_l(i) - v2 * ro_l(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      fro_l(i) = mx_l(i)
	   fmx_l(i) = vx * mx_l(i) + pr + b2 / 2.0d0 - bx_l(i) * bx_l(i)
	   fmy_l(i) = vx * my_l(i) - bx_l(i) * by_l(i)
	   fmz_l(i) = vx * mz_l(i) - bx_l(i) * bz_l(i)
	   fbx_l(i) = 0.0d0
	   fby_l(i) = vx * by_l(i) - bx_l(i) * vy
	   fbz_l(i) = vx * bz_l(i) - bx_l(i) * vz
	   fen_l(i) = (en_l(i) + pr + b2 / 2.0d0) * vx - (bx_l(i) * vx + by_l(i) * vy + bz_l(i) * vz) * bx_l(i)

      vx = mx_r(ip1) / ro_r(ip1)
      vy = my_r(ip1) / ro_r(ip1)
      vz = mz_r(ip1) / ro_r(ip1)
      b2 = bx_r(ip1) * bx_r(ip1) + by_r(ip1) * by_r(ip1) + bz_r(ip1) * bz_r(ip1)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_r(ip1) - v2 * ro_r(ip1) / 2.0d0 - b2 / 2.0d0) * gmm1
      fro_r(i) = mx_r(ip1)
	   fmx_r(i) = vx * mx_r(ip1) + pr + b2 / 2.0d0 - bx_r(ip1) * bx_r(ip1)
	   fmy_r(i) = vx * my_r(ip1) - bx_r(ip1) * by_r(ip1)
	   fmz_r(i) = vx * mz_r(ip1) - bx_r(ip1) * bz_r(ip1)
	   fbx_r(i) = 0.0d0
	   fby_r(i) = vx * by_r(ip1) - bx_r(ip1) * vy
	   fbz_r(i) = vx * bz_r(ip1) - bx_r(ip1) * vz
	   fen_r(i) = (en_r(ip1) + pr + b2 / 2.0d0) * vx - (bx_r(ip1) * vx + by_r(ip1) * vy + bz_r(ip1) * vz) * bx_r(ip1)

      ro_a(i) = (ro_l(i) + ro_r(ip1)) / 2.0d0
      mx_a(i) = (mx_l(i) + mx_r(ip1)) / 2.0d0
      my_a(i) = (my_l(i) + my_r(ip1)) / 2.0d0
      mz_a(i) = (mz_l(i) + mz_r(ip1)) / 2.0d0
      bx_a(i) = (bx_l(i) + bx_r(ip1)) / 2.0d0
      by_a(i) = (by_l(i) + by_r(ip1)) / 2.0d0
      bz_a(i) = (bz_l(i) + bz_r(ip1)) / 2.0d0
      en_a(i) = (en_l(i) + en_r(ip1)) / 2.0d0

      var_d(1, i) = ro_l(i) - ro_r(ip1)
      var_d(2, i) = mx_l(i) - mx_r(ip1)
      var_d(3, i) = my_l(i) - my_r(ip1)
      var_d(4, i) = mz_l(i) - mz_r(ip1)
      var_d(5, i) = bx_l(i) - bx_r(ip1)
      var_d(6, i) = by_l(i) - by_r(ip1)
      var_d(7, i) = bz_l(i) - bz_r(ip1)
      var_d(8, i) = en_l(i) - en_r(ip1)
   enddo

   call eigenspace (eigenvalue, eigenvector_l, eigenvector_r, ro_a, mx_a, my_a, mz_a, bx_a, by_a, bz_a, en_a, nx, gm)

   do i = 1, nx
      do i_dim = 1, 8
         var_c(i_dim) = sum (eigenvector_l(:, i_dim, i) * var_d(:, i))
      enddo

      var_c = var_c * abs (eigenvalue(:, i))

      fro(i) = (fro_l(i) + fro_r(i) + sum (var_c * eigenvector_r(:, 1, i))) * 0.5d0
      fmx(i) = (fmx_l(i) + fmx_r(i) + sum (var_c * eigenvector_r(:, 2, i))) * 0.5d0
      fmy(i) = (fmy_l(i) + fmy_r(i) + sum (var_c * eigenvector_r(:, 3, i))) * 0.5d0
      fmz(i) = (fmz_l(i) + fmz_r(i) + sum (var_c * eigenvector_r(:, 4, i))) * 0.5d0
      fbx(i) = (fbx_l(i) + fbx_r(i) + sum (var_c * eigenvector_r(:, 5, i))) * 0.5d0
      fby(i) = (fby_l(i) + fby_r(i) + sum (var_c * eigenvector_r(:, 6, i))) * 0.5d0
      fbz(i) = (fbz_l(i) + fbz_r(i) + sum (var_c * eigenvector_r(:, 7, i))) * 0.5d0
      fen(i) = (fen_l(i) + fen_r(i) + sum (var_c * eigenvector_r(:, 8, i))) * 0.5d0
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine roe
